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We calculate the superfluid stifltness of 2D lattice hard-core bosons at half-filling (equivalent to 
the S = 1/2 XY- model) using the squared winding number quantum Monte Carlo estimator. For 
X Ly lattices with aspect ratio L^/Ly = R, and Lx,Ly — > oo, we confirm the recent prediction 
[N. Prokof'ev and B.V. Svistunov, Phys. Rev. B 61, 11282 (2000)] that the finite-temperature 
stifi'ness parameters and determined from the winding number differ from each other and 
from the true superfiuid density ps . Formally, p^ ps in the limit in which Lx ^ oo first and then 
Ly ^ oo. In practice we find that p^ converges exponentially to p^ for R > 1. We also confirm that 
for 3D systems, p^ — p^ — pY = ps for any R. In addition, we determine the Kosterlitz-Thouless 
transition temperature to be Tkt / J = 0.34303(8) for the 2D model. 



I. INTRODUCTION 

In the usual Kosterlitz-Thouless^ description of a 2D 
superfluid,^ excitations of bound vortex-antivortex pairs 
renormalize the superfluid density Ps{T). Here, Ps{T) is 
the coefficient of the square of the gradient of the phase in 
the free energy expression that governs long-wavelength 
phase fluctuations. Above a critical temperature TKT^ 
the vortex pairs unbind and Ps{T) drops to zero. The 
well-known Nelson-Kosterlitz formula,'^ 



Tkt 



2 Ps {Tkt), 



(1) 



relates Tkt to the discontinuity in the superfluid den- 
sity Ps{Tkt) at Tkt- Using quantum Monte Carlo sim- 
ulations in a real-space basis, one can calculate a stiff- 
ness parameter p^ (T), related to Ps{T), using a squared 
winding number estimator.^ On an Lx x Ly torus, we can 
define two stiffness parameters (or helicity moduli) 



p^iT) ^ T^{W^), 



p^{T) = T-^{W^), 

where the integer winding numbers Wa, a 
defined according to 







(2) 
(3) 

x,y, are 
(4) 



where Jq, is the boson current operator and [3 = l/T. 
Typically, systems with aspect ratio R = L^/Ly = 1 
have been studied, and in the limit where Lx and Ly 
become infinite, one might expect that p^ {T) is equal 
to Ps{T). However, recently Prokof'ev and Svistunov^ 
have noted that in addition to the Kosterlitz-Thouless 
vortex-antivortex pairs, topological excitations present 
for an Lx x Ly torus lead to a renormalization of the 
stiffness parameters. In this case, at finite temperature. 



even in the limit Lx,Ly — > oo, p^ {T) and p^ (T) de- 
pend upon the aspect ratio and in general do not equal 
the superfluid density Ps(T).^ In particular, for Lx = Ly, 
pY {Tkt) / Ps{Tkt) ~ 0.9998. Hence, calculations of T/fr 
utilizing Eq. (1) and assuming p^ — ps can be expected 
to be affected by a small (typically negligible) systematic 
error. 

In this paper we present a quantum Monte Carlo study 
of the aspect ratio dependence of p^ and p^ , organized 
as follows. In Section II we review the argument of Ref. 5 
and summarize their conclusions. Section HI contains 
the results of our quantum Monte Carlo study on the 
S — 1/2 XY- model. In HI A we show the dependence 
of pY and p^ on R for 2D and 3D square lattices, and 
illustrate the need for an accurate, independent estimate 
of Tkt for this model to use as a benchmark. In Section 
IIIB we use the Weber-Minnhagen RG scaling relation^ 
on systems with R = 1 and 4 to obtain a very accurate 
estimate of Tkt- In Section IIIC we explore a method 
of calculating Tkt using the predicted aspect ratio de- 
pendence of I , and illustrate the need for precise 
finite-size scaling of the data in order to obtain satis- 
factory agreement with the benchmark Tkt- Once the 
finite-size scaling behavior of p^ j p^ is understood, it is 
straightforward to confirm the conclusions of Ref. 5 using 
our data. 




FIG. 1: A 2D superfluid with vortex excitations. The vortices 
are located at the points where the flux tube of circulation 
pierces the surface of the torus (-1- and — refer to a vortex 
and an ant i vortex). 
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II. STIFFNESS PARAMETERS 

Consider the Lx X Ly torus shown in Fig. 1. Imagine 
that its surface is coated with a 2D superfluid and that 
a tube of quantized circulation h/m penetrates the torus 
as shown in Fig. 1. There is an antivortex in the super- 
fluid layer at the point where this flux tube passes into 
the torus, and a vortex at the point where it leaves the 
surface. These excitations are the Kosterlitz-Thouless 
vortex-antivortex pairs. On a finite torus one can also 
envision a situation in which the flux tube moves through 
the cross-section and into the center of the torus as illus- 
trated in Fig. 2. In this case a single unit of quantized 
flux has entered the torus (Fig. 2(c)), and the excitation 
energy is 

2 



Ps{T) 



2tt 
Lx 



Lx L^, 



(5) 



Here, Ps{T) is the usual Kosterlitz-Thouless superfluid 
density which is renormalized by the vortex-antivortex 
pairs. Naturally, there arc excitations with energy e^;^^ 
associated with flux tubes containing i quanta. Similarly, 
there are excitations associated with a flux tube which 
threads around the inside of the torus along the Lx di- 
rection. The energies of these excitations will vary as 
Lx/ Ly. 




FIG. 2: Cross section of a torus, showing a flux line as it 
passes through to the center of the superfluid. In (a) the 
vortex-antivortex pair is maximally separated, in (b) the pair 
moves closer together and in (c) the pair has annihilated each 
other. The tube of quantized flux remains. 



At finite temperatures, the phase stiffness parameters 
and are affected by the vortex excitations. In the 
usual way, the stiffness is determined from the change 
in the free energy associated with an inflnitesimal flux (j) 
threading the torus. With (j) in the direction of the flux 
tube in Fig. 1, one obtains pY{T) from 



AFx = lp^iT)^J^ 



(6) 



with = 2TT(j)/{h/m). As discussed in Ref. 5, a calcula- 
tion of AF gives 

pY{T) = p,{T) ( 1 - hn^i-)x\ , (7) 



T Lx' 'V 
where £ is the number of quanta in the flux tube, and 



The stiffness in the y-direction is obtained by replacing 
Ly/Lx by Lx/Ly in Eqs. (7) and (8). 

If one lets Lx go towards infinity, keeping Ly finite, 



(Lx/Ly) T 
47rV,(T) ' 



(9) 



so that 



Px goes to zero as expected for a ID system at 
finite T. Alternatively, if the aspect ratio is such that Lx 
is large compared to Ly, then 



Ps{T) 



Sn'^PsiT) Lx -2-n- 
T L,. 



2 Ps(T) 



(10) 

and the difference between {T) and Ps{T) vanishes 
exponentially as Lx/Ly increases. Finally, in 3D the ad- 
ditional factor of Lz which occurs in Eq. (7) assures the 
convergence of p^ to ps for all R as the size of the system 
goes to infinity. 

In the 2D superfluid, the jump in the stiffnesses at the 
Kosterlitz-Thouless temperature also clearly depends on 
the aspect ratio. At Tkt, using the Nelson-Kosterlitz 
relation (1), one has'"' 



{Tkt) = Ps{Tkt) (l - ^ , (11) 



with 



(12) 



In a similar way, p^ {Tkt) and {(.'^)y are obtained by 
replacing Ly/Lx by Lx/Ly in the above two equations. 
Hence, one can determine the ratios / ps and p^ / Ps 
at Tkt by evaluating the sums for and Ta- 
ble I shows the result of this calculation for various aspect 
ratios. 



TABLE I: Stifi'ness parameter ratios at Tkt, for different as- 
pect ratios R = Lx/Ly. The entry * was evaluated exphcitly,* 
while the remaining entries were evaluated to arbitrary nu- 
merical accuracy using Waterloo's Maple V software. 



R 


W 1 
Px /ps 




W 1 
Py /Ps 




1 


0.9998247 




0.9998247 




2 


0.9532407 


1 - 


1.222613 X 10- 


-9 


3 


0.7533903 


1- 


6.395505 X 10" 


15 


4 


1/2* 


1 - 


2.973776 x 10" 


20 


5 


0.2977676 


1 - 


1.296322 X 10" 


2B 


6 


0.1663429 


1 - 


5.424861 X 10" 


31 


7 


0.0893401 


1 - 


2.207141 X 10" 


36 


8 


0.0467593 


1 - 


8.796634 x 10" 


42 
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III. MONTE CARLO CALCULATIONS 

In order to study this aspect ratio dependence of 
and we employ a model of a two dimensional super- 
fluid using a hard-core boson Hamiltonian at half filling, 
which is equivalent to the quantum 5=1/2 XY model 
defined by 

H=^jY, [S^S^ + SfS^^] . (13) 

Here Sf and Sf are the x and y components of a spin 

1/2 operator at site i, and the sum is carried out over all 
nearest- neighbor spin pairs This model is known 

to have a Kosterlitz-Thouless type transition from pre- 
vious Monte Carlo simulations carried out on square 
lattices. We carry out simulations using the stochastic 
scries expansion (SSE) quantum Monte Carlo method, -'^'^ 
that has previously been applied to this and other spin 
and boson models. The basis of the SSE method is im- 
portance sampling of the power series expansion of the 
partition function: 

Z = Tr{e-^"} = J2Y1 ^^m-b), (14) 

7 n=0 

where (3=1 /T is the inverse temperature and the trace 
has been written as a sum over diagonal matrix elements 
in a basis {I7)}. In our case this is the standard basis of 
z spin quantization. In the simulations we employ an ef- 
ficient cluster-type method for sampling the terms of the 
expansion, called the directed-loop algorithm,^^ which in 
the present case of zero field is similar to loop algorithms 
previously employed on the spin 1/2 XY model. Note 
that all contributing expansion-orders n in (14) are sam- 
pled, and all results are exact within statistical errors. 

A direct estimator for the superfluid (spin) stiffnesses is 
identified in the SSE quantum Monte Carlo in a way anal- 
ogous to world-line quantum Monte Carlo methods. 
Starting from the definition Eq. (6) and taking the sec- 
ond derivative of the free energy (per spin) with re- 
spect to a twist 9 in the boundary condition will lead 
to the expressions given in Eq. (2) and (3). The wind- 
ing number in these equations is now a measure of the 
net spin currents flowing around the periodic system, 
Wa = {N^ - N^)/La, {a = x,y) where N^'^ is the 
number of operators transporting spin to the "right" or 
"left" along the a direction in the SSE configuration. 



A. Spin Stiffness in 2D and 3D 

In order to test the predictions of Prokof'ev and 
Svistunov,^ we carried out a series of simulations using 
the Hamiltonian, Eq. (13), for systems with various as- 
pect ratios in two and three dimensions. In general our 
calculations confirm the predictions of Ref. 5, namely 




T/J 

FIG. 3: The spin stiffness order parameter Eqs. (2), (3) cal- 
culated for a (a) 16x16x64 system and a (b) 64x16 system. 
The long lattice direction is z and x, respectively. The slight 
deviation of pY (T) around T/J « 1.0 in (a) disappears in 
the limit of largo system size. The thick dashed line in (b) 
represents the equation = 2T/nJ, and the thin vertical 
dotted line marks the corresponding estimate of Tkt/J- 



that in the 3D system, the difference between the com- 
ponents of p'^ {T) {a = x,y,z) vanishes exponentially, 
while in 2D, p^{T) can differ significantly from p^(T) 
depending on the value of the aspect ratio R. Specific ex- 
amples are illustrated in Fig. 3. As we see in Fig. 3(a), in 
the 3D case, for LjjLy = 1 and L^/ Ly = A the three dif- 
ferent spin stiffness parameters converge to one value. In 
contrast, for the 2D case (Fig. 3(b)), with R = L^/Ly = 4 
the spin stiffness in the long direction is significantly less 
that the spin stiffness in the short direction for a large 
range of T. As we will see below, in 2D p^{T) does not 
converge to (T) in this temperature range even in the 
limit of large system sizes (over 6.5 x 10** spins). 

A rough estimate of Tkt and the ratio 
pf {Tkt) / {Tkt) can be obtained from Fig. 3(b), by 
drawing a straight line p'^ = 2T/7rJ. The point where 
this line intersects w ps is a finite-size measurement 
of Tkt/J, and the ratio p^ / p^ at this T can be 
compared to the value 1/2 from Table I (which we 
expect to be exact in the limit of infinite system size). 
It is clear from this simple demonstration that if one 
wishes to study the spin stiffness dependence on R in 
detail, an accurate independent estimate of Tkt will 
facilitate quantitative comparisons of p^ j p'y with the 
entries in Table I. The next section of this paper is 
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therefore dedicated to the measurement of Tkt- 



B. Weber-Minnhagen'^ determination of Tkt 

Wc now give details on the independent estimate for 
Tkt which will be used as a benchmark for Txt{L — > oo) 
in the discussion to follow. Our approach follows that 
of Harada and Kawashima,^° who carried out simula- 
tions of the quantum S = 1/2 XY model on nine square 
lattices of size L = 8 to 128. Their analysis found a 
Tkt = 0.3427(2) J, based on a Webcr-Minnhagen^ scal- 
ing fit for the system size dependence of the spin stiffness. 



Specifically, this scaling form states that as T — > 



KT' 



Ps{L) = Ps(oo) 



1 + 



1 



1 



2 ln(L) + C) 



(15) 



only, should have a minimum in s-t Tkt/ J- This is 
illustrated in Fig. 4. In general, we find that the qual- 
ities of the curves can differ significantly depending 
on the "lattice sequence" included in the fit. An accept- 
able lattice sequence should be one that produces a 
minimum of approximately unity. We therefore system- 
atically eliminate the smallest Ly values from the lattice 
sequence until Xmin ~ 1- With our tempering data, this 
occurs when the largest excluded data point is Ly = 8, 
however to be cautious we also excluded Ly = 12. Quali- 
tatively, the shift in xLin is very small (AT/ J < 0.0001) 
when Ly > 16 points are excluded from the data set, 
which in part is due to the trivial statistical systematic 
shift resulting from the elimination of data points from 
the fit. Finally, the curves have little quantitative de- 
pendence on the maximum Ly included, as long as this 
Ly > 64. 



where PsiL) now represents the spin stiffness at a finite 
lattice with a linear size of L spins, Ps{oo) is given by 
Eq. (1), and C is some unknown constant. One can also 
reformulate Eq. (15) in the linear form. 



Y{L) =Trln{L) + C, 



where 



Y{L) = 



P^jL) 1 2 
F{R) T tt 



(16) 



(17) 



and C has been re-scaled to absorb some simple con- 
stants. The factor F{R) has been included in order to 
correct according to ps = p"^ /F{R) for square sys- 
tems, where F{1) = 0.9998247 from Table I. 




0.342 



0.344 



FIG. 4: values (per degree of freedom) for R = 1, lattice 
sequence L = 16, 24, 32, 40, 48, 56, 64, 72, 96, 128 and i^ = 4, 
Ly = 16, 20, 24, 32, 40, 48, 56, 64. 



In order to obtain results on a dense temperature grid, 
we did an extensive series of simulated tempering^'^ sim- 
ulations within a T range centered about the approx- 
imate Tkt ~ 0.343J.i° Fits to Eq. (16), adjusting C 
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FIG. 5: Simulated tempering data for py , where each data 
point represents 2 x 10^ Monte Carlo production steps. Statis- 
tical errors are at most of the order of the symbol sizes. The 
star on the vertical axis is {Ly — » oo) determined from 
Tkt ~ 0.343J using Eq. (1). Inset: data scaled to the 
form Eq. (16). Data included in the fit is plotted with solid 
symbols, and data for smaller system sizes that was omitted 
from the x^ At is plotted with open symbols. The solid lines 
have constant slopes of tt. 



Using the scaling form Eq. (16) with the data illus- 
trated in Fig. 5, we obtain three estimates for Tkt(oo): 
for (i) i? = 4, using F{R) = 1, (ii) i? = 1, with the ap- 
proximation F(R) = 1, and (iii) R = 1, with the correct 
factor F{R) = 0.9998247 from Table I. The correspond- 
ing results for the transition temperatures are t'^^j,/ J = 

-4 rp{ii 



0.34302 ± 1.4 X 10 



j^^/J = 0.34302 ± 0.9 X 10"* 

/J = 0.34304 ± 1.0 X 10-*, respectively, where 
quoted errors are one standard deviation. From this, we 



and T, 



(Hi) 



KT 



see that T, 



ii) 



KT 



T 



(■«) 
KT 



Tkt^ to within our statistical er- 
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ror. To obtain a more accurate estimate for the transition 
temperature, we may also perform a weighted average of 
our vahics for Tl^j, and TjjV, to get Tkt/J = 0.34303(8). 
This value agrees to the previous^° most accurate Tkt = 
0.3427(2) J to within two standard deviations. 

The Monte Carlo data and the fit at our estimated 
Tkt are shown in Fig. 5. As evident there, the unsealed 
data at T/J = 0.34303 do not clearly approach the 
infinite size limit determined with Eq. (1). However, 
when properly scaled according to Eq. (17), Y{L) does 
become linear according to the expected form Eq. (16) 
(Fig. 5 inset). 



C. The Spin Stiffness Ratio p^/pY at Tkt 

In this section we study the aspect ratio dependence 
of the spin stiffness in more general terms, and outline 
a procedure for calculating Tkt using p^ j and the 
results of Table I. Our aim here is not to obtain a more 
precise value for Tkt^ but to test how closely the pre- 
dicted values for the ratios can be observed in practice. 

To begin, we collected detailed data for p}^ (T) and 
p^iT) for 20 to 24 values of T/J between 0.325 and 
0.375. Simulations produced between 1 x 10^ and 5 x 10^ 
Monte Carlo steps at each T (depending on system size 
and statistical error bars), except for the very largest sys- 
tems which had an order of magnitude less Monte Carlo 
steps. System aspect ratios included were R =1,2 and 
4. Continuous data over T was obtained by fitting third- 
order polynomial curves. 

Fig. 6 shows the spin stiffness ratios p^ j calculated 
with this data, at Tkt I J = 0.34303, as was determined 
independently in Section IIIB above. As evident in the 
figure, this type of simple analysis is not very satisfy- 
ing, as the spin stiffness ratios do not obviously approach 
the values expected from Table I. Clearly there are sig- 
nificant finite-size effects at play, as seen previously in 
the Weber-Minnhagen scaling analysis of Section IIIB 
(Fig. 5). One way to quantitatively account for the scal- 
ing effects in our p^ j will now be outlined. To do 
so, we define a new measurement technique which can 
be used to estimate the Kosterlitz-Thoulcss transition 
temperature: namely, Tkt is defined as the temperature 
at which the ratio p^ / ps (in practice p^ / p^), passes 
through the corresponding value in column 2 of Table I 
(for R —2 and 4 only). Additionally, for the purposes 
of comparison, Tkt can be obtained by calculating the 
crossing of Eq. (1) with the p^ (T) data for each system 
size Lx X Ly. The results of these two procedures are 
summarized in Fig. 7. 

Several interesting observations can be made from this 
figure. First, if one were only to consider smaller system 
sizes {Ly < 64), the data would appear to consistently 
approach an infinite size limit of Tkt/J ~ 0.35. How- 
ever, at Ly K, 64 we observe that all of the data sets in 
Fig. 7 undergo subtle changes in slope that could indeed 
suggest a significantly lower value for Tkt in the infinite 
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FIG. 6: The spin stiffness ratios for R = 2 and i? = 4 at 
T/J = 0.34303. At T = Tkt and L:,,Ly oo, the values 
are expected to approach / p^ = 0.95324 and 1/2, respec- 
tively. 
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FIG. 7: The Kosterlitz-Thouless transition temperature cal- 
culated with the Nelson-Kosterlitz formula (NK, Eq. (1)), and 
with the spin stiffness aspect ratios (PS, Table I), as a func- 
tion of inverse linear system size. The infinite size limit is 
Tkt ~ 0.343 J (see Section lllB of text) and is plotted as a 
star on the vertical axis. 



size limit. Note in particular the change in the sign of 

the slope for the PS, R — 2 fine at Ly « 64. Without 
any a priori knowledge of the form of the scaling laws 
for the curves in Fig. 7, it would be difficult to conclude 
whether they approach the Kosterlitz-Thoidcss tempera- 
ture of Tkt ~ 0.343 J calculated in Section lllB, and in 
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Ref. 10, illustrated on the vertical axis of Fig. 7. 
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FIG. 8: The Kosterlitz-Thoulcss transition temperatures 
scaled according to Eq. (19). In (a), the Tkt data for 7? = 4 
becomes approximately linear (excluding the smallest system 
size, Ly = 8) for Lo « 0.01. In (b), all Tkt curves become 
approximately linear for Lo in a narrow range around 0.3. 
Straight lines are included as guides for the eye. The infinite 
size limit is Tkt = 0.343 J (see Section IIIB of text) and is 
plotted as a star on the vertical axis. 



We now consider a scaling law for Tkt that arises from 
the Kosterlitz RG scaling equations. ^'^ As discussed in 
Ref. 13, the functional form of the correlation length sug- 
gests that, as T — > Tkt, £, diverges according to 



(18) 



where t = {T — Tkt) /Tkt > and 6 is a constant. 
If we identify ^ = L/Lq {Lq some microscopic length), 
and T = Tkt{L), then we can write the system size- 
dependent Kosterlitz-Thouless transition temperature in 
terms oi L/Lq as 



Tkt{L) = Tkt{oo) 



\n\L/Lo) 



(19) 



Fig. 8 illustrates the Tkt(L) data of Fig. 7, rc-scalcd 
to the form of Eq. (19), for some value of Lq which ap- 
proximately linearizes the data. As this figure shows, 
by choosing appropriate values for Lq it is possible to 
obtain linear fits to Eq. (19) for all data sets save one. 
The exception is the spin stiffness aspect ratio crossing 
for i? = 2 (PS, i? = 2 in Figs. 7 and 8). In this case. 



the change in slope of the data set at Ly « 64 precludes 
any direct fitting to the scaled form Eq. (19). However, 
the data for Ly > 64 show rough evidence of this scaling 
trend, and it is therefore very likely that additional data 
for Ly > 128 would approach this Tkt{oo) consistently. 

As Fig. 8 shows, the approximate intercept of these 
scaled data sets are in good agreement with the value 
of Tkt{oo) calculated in Section IIIB using the Weber- 
Minnhagcn^ scaling form for the spin stiffness. In par- 
ticular, the consistency of the intercept for the data set 
PS, i? = 4, adds confidence to the assertion that Tkt{L) 
calculated using the spin stiffness aspect ratios is a good 
estimator for the Kosterlitz-Thouless transition temper- 
ature, consistent with Tkt calculated using other meth- 
ods. Thus the entries in Table I indeed appear to be accu- 
rate. Clearly, however, our value of Tkt/ J = 0.34303(8) 
obtained from the Webcr-Minnhagen scaling is the most 
accurate because it is based on a one-parameter fit, 
Eq. (16), as opposed to the two-parameter fit, Eq. (19). 



IV. CONCLUSIONS 

In conclusion, we have confirmed the prediction of 

Prokof'ev and Svistunov^ that in 2D the squared 
winding-number estimates of the finite-temperature spin 
stiffnesses and differ from each other for i? 7^ 1, 
while in 3D systems = = pY for ^-^ly R- 2D 
we also find that p^ approaches the Nelson-Kosterlitz su- 
perfluid density ps exponentially for i? > 1 as Lx,Ly — > 
00, whereas pjf approaches zero as i? — > 00 exactly as 
predicted. 

As illustrated in Fig. 8, the Tkt predicted from the 



ratios p^ / plj approach the Tkt predicted using other 
methods, substantiating the validity of the results listed 
in Table I. Neglecting the small differences between ps 
and p'^y on 2D 1/ X L (i? = 1) lattices is insignificant 
compared to presently typical statistical uncertainties. 
However, with any further increase in accuracy relative 
to that achieved in the data presented here, the effect has 
to be taken into account in order to avoid a systematic 
error. 
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